module frozen_flow
   !
   ! SUBROUTINES:
   !
   ! 1) frozen_cp_gamma
   ! 2) frozen_Mach
   ! 3) frozen_mi_k
   ! 4) frozen_inlet_pressure
   !
   ! Last update: 23 Jan 2012

   use user
   use Gibbs_1p4

contains

   !****************************************************************************

   ! Subroutine 01

   subroutine frozen_cp_gamma (nx, ny, T, Yi, cmod, itimax, itemax, tol_ej, &
         tol_n, po, OF, cp, gcp) ! The last two are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real + ghost)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real + ghost)
      integer, intent(in) :: cmod   ! Chemical model
      integer, intent(in) :: itimax ! Maximum number of iterations to solve e(j) - dissociation rate for j-reaction
      real(8), intent(in) :: tol_ej ! Tolerance for e(j) evaluation
      integer, intent(in) :: itemax ! Maximum number of iterations to solve n - the number of mols of products
      real(8), intent(in) :: tol_n  ! Tolerance for n evaluation
      real(8), intent(in) :: po     ! Stagnation pressure in the chamber
      real(8), intent(in) :: OF     ! Oxidant/Fuel ratio
      real(8), dimension(:), intent(in) :: T     ! Temperature (K)
      real(8), dimension(:,:), intent(in) :: Yi    ! Mass fraction (non-dim.)
      real(8), dimension(:), intent(out) :: cp    ! Specific heat at const pressure (J/kg.K)
      real(8), dimension(:), intent(out) :: gcp   ! gcp = gamma = Cp/Cv at the center of CV P (non-dim.)

      integer :: i

      do i = 1, nx*ny

         call GIBBS_EQUILIBRIO_dados_recebe (1, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T(i), po/101325.0d0, OF)
         call GIBBS_MACH1D_recebe1 (Yi(i, :))

         call GIBBS_cp_calculo

         call GIBBS_gama_calculo

         call GIBBS_MACH1D_envia1 (cp(i), gcp(i))

      end do

   end subroutine frozen_cp_gamma

   !****************************************************************************

   ! Subroutine 02

   subroutine frozen_Mach (Rgp, gcp, u, T, M) ! The last one is output
      implicit none
      real(8), dimension(:), intent(in) :: Rgp  ! Gas constant (J/kg.K)
      real(8), dimension(:), intent(in) :: gcp  ! gcp = gamma = Cp/Cv at the center of CV P (non-dim.)
      real(8), dimension(:), intent(in) :: u    ! velocity (m/s)
      real(8), dimension(:), intent(in) :: T    ! temperature (K)
      real(8), dimension(:), intent(out) :: M    ! Mach number

      M = u/dsqrt(gcp*Rgp*T)

   end subroutine frozen_Mach

   !****************************************************************************

   ! Subroutine 03

   subroutine frozen_mi_k(nx, ny, T, v_type, visc_cte, k_type, k_cte, vlp, kp) ! The last two are output
      implicit none
      integer, intent(in) :: nx ! Number of volumes in csi direction (real + ghost)
      integer, intent(in) :: ny ! Number of volumes in eta direction (real + ghost)
      integer, intent(in) :: v_type   ! v_type = 1 -> constant; v_type = 2 -> frozen
      real(8), intent(in) :: visc_cte ! constant dynamic viscosity (only for v_type = 1)
      integer, intent(in) :: k_type   ! k_type = 1 -> constant; k_type = 2 -> frozen
      real(8), intent(in) :: k_cte    ! constant thermal conductivity (only for k_type = 1)
      real(8), dimension(:), intent(in) :: T   ! temperature (K)
      real(8), dimension(:), intent(out) :: vlp ! dynamic viscosity (Pa.s)
      real(8), dimension(:), intent(out) :: kp  ! thermal conductivity (W/m.K)

      integer :: i

      if (v_type == 2 .or. k_type == 2) then

         do i = 1, nx*ny

            call GIBBS_MACH1D_recebe2 (T(i))

            call GIBBS_CONGELADO_calculos2

            call GIBBS_CONGELADO_envia2 (vlp(i), kp(i))

         end do

      end if

      if (v_type == 1) vlp = visc_cte

      if (k_type == 1) kp = k_cte

   end subroutine frozen_mi_k

   !****************************************************************************

   ! Subroutine 04

   subroutine frozen_inlet_pressure (gamma, Rg, p_ch, u_in, T_in, p_in)

      real*8, intent(in) :: gamma    ! gamma = Cp/Cv at the center of CV P (non-dim)
      real*8, intent(in) :: Rg       ! Gas constant (J/kg.K)
      real*8, intent(in) :: p_ch    ! pressure at combustion chamber (Pa)
      real*8, intent(in) :: u_in     ! inlet velocity (m/s)
      real*8, intent(in) :: T_in     ! inlet temperature (K)

      real*8, intent(out) :: p_in  ! inlet pressure (Pa)

      real*8 :: k    ! auxiliar
      real*8 :: M_in ! inlet Mach number

      M_in = u_in/dsqrt (gamma*Rg*T_in)

      k = 1.0d0 + 0.5d0*(gamma - 1.0d0)*(M_in**2)

      p_in = p_ch/(k**(gamma/(gamma - 1.0d0)))

   end subroutine frozen_inlet_pressure

   !****************************************************************************

end module frozen_flow
